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EXECUTIVE SUMMARY 


No closed form solutions exist for the elastic-plastic /-integral for surface cracks due to the 
nonlinear three-dimensional (3D) nature of the problem. Traditionally, each surface crack case must 
be analyzed with a unique and time-consuming nonlinear finite element analysis (FEA). Addition- 
ally, knowledge of nonlinear FEA, plasticity theory, and elastic-plastic fracture mechanics is required 
to reliably execute this type of analysis. To simplify this process, the authors have developed and ana- 
lyzed an array of 600 3D nonlinear finite element models for surface cracks in flat plates under ten- 
sion loading. The solution space covers a wide range of crack shapes and depths (shape: 0.2 <a/c< 1, 
depth: 0.2<u/i?<0.8) and material flow properties (elastic modulus to yield ratio: 100 < El o vs < 1,000, 
and hardening: 3 < « < 20). The solution of this large array of nonlinear models was made practical 
by computer routines that automate the process of building the finite element models, running the 
nonlinear analyses, post-processing model results, and compiling and organizing the solution results 
into multidimensional arrays. The authors have developed a methodology for interpolating between 
the geometric and material property variables that allows the user to estimate the /-integral solution 
around the surface crack perimeter as a function of loading condition from the linear-elastic regime 
through the elastic-plastic regime. In addition to the /-integral solution, the complete force versus 
crack mouth opening displacement ( CMOD ) record is estimated. The solution space and interpo- 
lation routines have been extensively verified throughout the range of applicable geometries and 
materials. The user of this interpolated solution space need only know the crack and plate geometry 
and the basic material flow properties to reliably evaluate the full surface crack /-integral and force 
versus CMOD solution; thus, a solution can be obtained very rapidly by users without elastic-plastic 
fracture mechanics modeling experience. Enabling a convenient surface crack solution in the fully 
elastic-plastic domain, which up to now was only available for linear-elastic conditions, will signifi- 
cantly reduce the costs typically associated with evaluating elastic-plastic surface crack behavior. The 
authors hope this will promote advanced surface crack testing methodologies in the elastic-plastic 
regime and simplify aspects of surface crack structural assessment using the elastic-plastic /-integral. 
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ELASTIC-PLASTIC J-INTEGRAL SOLUTIONS LOR SURFACE CRACKS IN TENSION 
USING AN INTERPOLATION METHODOLOGY 

1. INTRODUCTION 


Surface cracks are among the most common defects found in structural components and 
frequently reach failure once the limits of linear-elastic fracture mechanics (LEFM) have been 
exceeded and elastic-plastic fracture mechanics (EPFM) governs. For example, in the aerospace 
industry EPFM conditions commonly exist due to small critical crack sizes, thin walls, and high 
loading conditions inherent to lightweight, high-performance structures. Conversely, in the nuclear 
and petroleum industries EPFM conditions often prevail due to the use of lower strength, very high- 
toughness materials. EPFM assessment has become significantly more accessible through improved 
Unite element interfaces such as FEACrack™ or ABAQUS® CAE, 1 ’ 2 but unfortunately the cost of 
such assessments in analysis time remains a significant impediment to common use. Similar difficul- 
ties arise in assessing laboratory fracture toughness tests with surface cracks. In these tests, due to 
practical specimen size limitations, the material fracture toughness is commonly not reached until 
well beyond the LEFM limit. Surface crack fracture testing is hindered significantly by the lack of 
a readily available set of surface crack solutions for the nonspecialist to evaluate the elastic-plastic 
/-integral, crack mouth opening displacement ( CMOD ) values, or deformation state of a test speci- 
men at fracture. A convenient set of elastic-plastic surface crack solutions could help mitigate many 
of these obstacles. 

Figure 1 illustrates a typical semielliptical surface crack in a flat plate in tension, where W is 
the total width of the plate, L is the length, B is the thickness, 2c is the total surface crack length, 
and a is the crack depth. The parametric crack front angle, <j>, is defined in the conventional form as 
shown in section A-A of figure 1 where 0 = 0° and 180° at the free surface, and (j) = 90° at the deepest 
point of the crack. The CMOD is defined as the axial displacement between the two crack mouth 
faces on the free surface of the plate at the centerline of the crack. A uniform far-field displacement, 
8, or stress, o, is designated by the arrows on the top and bottom faces of the plate. 


r- CMOD 



|-« c ► 

2c ► 

Section A-A 


Figure 1. Semielliptical surface crack in a flat plate. 

The Newman-Raju (NR) equations 3 ’ 4 provide a robust and convenient linear-elastic stress 
intensity ( Kj ) solution for surface cracks in tension and bending. Other common linear-elastic codes, 
such as NASGRO® utilize an interpolated table look-up approach for surface crack X / solutions. 5 
Unfortunately, no concise solutions exist for elastic-plastic surface cracks due to their nonlinear 
three-dimensional (3D) nature. The contemporary method for solving the elastic-plastic surface 
crack problem requires a unique and time-consuming nonlinear finite element analysis (FEA) for 
each surface crack geometry of interest. To reliably execute this type of analysis requires knowl- 
edge of nonlinear FEA, plasticity theory, and EPFM. In the years prior to the advent of routine 
finite element based fracture mechanics analysis, many researchers provided alternative and robust 
engineering solutions to the elastic-plastic surface crack problem, though subject to many practical 
limitations. An excellent summary of the development of elastic-plastic /-integral solutions up to the 
year 1999 is given by McClung et al. 6 Apart from FEA, the commonly used methods for calculating 
elastic-plastic /-integral solutions usually follow one of two basic techniques: 7 the Electric Power 
Research Institute (EPRI) approach 8-10 or the reference stress method (RSM). 11 The EPRI method 
assumes a Ramberg-Osgood 12 material stress-strain behavior, requires an estimate for the yield load, 
and requires the use of nondimensional h 1 influence factors in the calculation of the /-integral. The 
h] influence factors are a function of both constitutive relationship and geometry and are obtained 
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from a fully plastic elastic-plastic finite element solution. Tables of /q factors for surface cracks in flat 
plates have been generated by several researchers 6 ’ 13-1 5 for limited combinations of material prop- 
erties, geometry values, and surface crack front angle, <p. Yagawa et al. 13 also developed influence 
factors for the fully plastic crack opening displacement for the geometries and materials covered in 
their study. In addition, Kim et al. 16 performed investigations into the RSM for a limited number of 
surface crack geometries and Ramberg-Osgood material combinations. 

The EPRI and RSM techniques have found wide application in analysis of structures but 
have limited application in the assessment of surface crack laboratory tests. To understand with 
sufficient precision the crack front conditions at the point the fracture toughness is reached in 
an experimental surface crack test requires knowledge of the force, P, versus CMOD response, the 
elastic-plastic deformation state of the specimen, and detailed knowledge of the /-integral versus (p 
relationship as a function of deformation. The current RSM and EPRI solutions for surface cracks 
do not provide the user with the full P versus CMOD trace which serves as the most fundamental 
connection between experiment and analysis. The CMOD value provides the most robust predictor 
of the /-integral values at the crack tip. 7 ’ 17 Most conventional fracture toughness test specimens 
such as the compact tension, C(T), or single edge notched bend are considered two-dimensional (2D) 
in nature and report a single value of fracture toughness representing the average driving force along 
the entire crack front. Conversely, the surface crack toughness test is highly 3D, and the toughness 
values are usually reported as a single, local value of toughness at a given (p location along the crack 
front. Most of the current RSM and EPRI solutions only have solution values at a limited number of 
(p locations. In addition, the EPRI and RSM techniques have / versus (p relationships that are based 
on either linear-elastic solutions (RSM) or fully plastic solutions (EPRI) which do not capture the 
changes in the / versus (p distribution and maximum /-integral location with increasing deformation. 

To overcome these shortcomings and provide a simple and robust method for analyzing sur- 
face crack tension tests, the authors have developed and analyzed an array of 600, 3D nonlinear 
finite element models (FEMs) for surface cracks in flat plates under tension loading. The solution 
space covers a wide range of crack geometric parameters and material properties. The solution of 
this large array of nonlinear models was made practical by computer routines that automate the pro- 
cess of building the FEMs, running the nonlinear analyses, post-processing model results, and com- 
piling and organizing the solution results into multidimensional arrays. The authors have developed 
a methodology for interpolating between the geometric and material property variables that allows 
the user to estimate the /-integral solution around the surface crack perimeter (< p ) as a function of 
loading condition from the linear-elastic regime through the elastic-plastic regime. In addition to the 
/-integral solution, the complete force versus CMOD record is estimated. The user of this interpo- 
lated solution space need only know the crack and plate geometry and the basic material flow prop- 
erties to reliably evaluate the full surface crack /-integral and force versus CMOD solution; thus, 
a solution can be obtained very rapidly by users without EPFM modeling experience. 
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2. COMPUTATIONAL PROCEDURES 


This project is logistically intense. Though computationally, each part of the process of build- 
ing this new solution space follows mostly well-established paths, combining those parts effectively 
into a functional whole requires planning at every level. This section provides insight into the basic 
computational procedures used in the FEMs which constitutes the solution space as well as the 
parameters of the geometric and material property variables that define the overall solution space. 
Briefly discussed in the conclusion are the logistics of building, executing, and then assembling the 
solution space — made practical only through automation. 


The constitutive model for this study employs incremental von Mises plasticity theory in 
a conventional small geometry change (small strain) setting. The material is assumed to be isotronic 
and homogenous with flow behavior following a linear plus power law (LPPL) representation of 
the stress-strain response. For the LPPL model, the uniaxial stress-strain curve follows a linear then 
power law model given by 


where e is engineering strain, n is the strain hardening exponent, a is engineering stress, a is a rep- 
resentative yield stress, and e ys is a corresponding yield strain defined by £ ys = o ys /E, with the elastic 
modulus, E. For this analytical work, o ys is equivalent to the proportional limit and defines the lim- 
its of the initial linear-elastic portion of the response. This approach to modeling the stress-strain 
response has numerous advantages in this study: first, the model lends itself to easy implementation 
in the planned material space because the curve is fully defined by three parameters (o ys , e ys or E, and 
n ); second, the model faithfully represents the stress-strain response of a large majority of structural 
metals; and third, unlike the commonly used Ramberg-Osgood model, prior to exceeding the yield 
strain, the LPPL model is purely linear-elastic, and thus does not prematurely accumulate plastic 
strain. Unlike the research FEA code utilized for this study, WARP3D, 18 most commercial finite ele- 
ment codes do not natively support the LPPL model for plasticity analysis. Therefore, to make the 
analysis methods used in this study more universally repeatable in other analysis codes, the smooth 
LPPL representation of the material response for £>£ ys was discretized into 20 line segments for use 
in the commonly available, multilinear stress-strain curve representation. In this way, the common 
and convenient way of describing the plastic hardening behavior through the hardening exponent, 
n, could be retained but the analysis performed in a way that is not dependent on the availability of 
a native LPPL model. The plastic strain, £ p[ , increments for each multilinear line segment were gra- 
dated from initially fine increments of 0.001 for eight line segments, to 0.005 for four line segments 
and then to 0.01 for the remaining seven line segments to ensure a smooth transition from linear- 
elastic to nonlinear behavior and to smoothly capture the power law hardening behavior. Figure 2 
illustrates the multilinear representation of the LPPL plastic response for a material with a yield 
strain, e ys = 0.003 (normalized elastic modulus of £70" =300) and various n values. 


2.1 Constitutive Model 
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Figure 2. Multilinear representation of the LPPL material for E!a ys = 300 
and various values of n. 


2.2 Solution Space 

The solution space for this array of models is four-dimensional. Two dimensions are used 
to describe surface crack geometric variation, and two dimensions are used to describe material 
property variation. The inclusion of nonlinear materials uniquely separates this work from previ- 
ous surface crack solutions for Kj, the material independent linear-elastic stress intensity factor. The 
material and geometric spaces must be carefully crafted to provide sufficient coverage for most com- 
mon engineering problems without becoming so large as to be intractable. The following sections 
discuss the choices and reasoning for the material and geometric dimensions of the solution space. 

2.2.1 Material Space 

As discussed above, the LPPL model for the stress-strain curve is fully defined by just three 
parameters. If the yield strength is normalized to unity for all materials ( cr v . v = 1), then only e and n 
are required to define the shape of the stress-strain curve throughout the space. For convenience of 
eliminating small fractional numbers, the reciprocal of the yield strain is commonly used, Elo ys . 
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Figure 3 illustrates the material space for the study described in terms of the six Elo ys and 
five n values resulting in 30 different material combinations. In all cases, o ys = 1 and Poisson’s ratio, 
v = 0.3. The names of several common engineering materials are overlaid on the material matrix in 
figure 3 to illustrate how some common materials are represented in the material matrix. The low 
E/(7 vs values of 100 to 200 are materials capable of high values of elastic strain, thus they have low 
elastic modulus and relatively high yield strength, such as many high-performance titanium and alu- 
minum alloys. The opposite end of the Elo space with values of E/o y = 1 ,000 have very little elastic 
strain capability due to high elastic modulus and low yield strength. Austenitic stainless steels are 
a common example of this material class. To practitioners not accustomed to considering materials 
based on the E/o vs ratio, it is important to realize that, for a fixed hardening exponent, n, the nonlin- 
ear material response is governed by this ratio. Materials that may not seem alike from an engineering 
perspective are actually equivalent in their nonlinear response. For example, consider three materials 
with similar n and equivalent yield strains such that El o ys = 400 in each case. This could easily be 
an aluminum alloy with yield strength around 180 MPa, a titanium alloy with yield strength around 
275 MPa, or a steel alloy with yield strength around 500 MPa. The normalized nonlinear response is 
the same in each case. This will be discussed further in a following section on normalization schemes. 
The effect of E/a ys on the load, P, versus CMOD behavior (and the resulting /-integral values) is 
strongest for lower values of Ela vs \ therefore, the authors chose smaller increments of Elo for 
Elo < 300 to provide more uniform coverage over the solution space. 
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Figure 3. Material space. 


The other dimension of the material space is the strain hardening exponent, n. The values of 
n range from 3 to 20, spanning the hardening characteristics of most all structural metals from very 
high strain hardening (n = 3) to almost elasticly perfect plastic behavior (n = 20). The specific values 
of n for this study were chosen to uniformly divide the strain hardening response in the a versus e p j 
space as shown in figure 2. 

2.2.2 Geometric Space 

Figure 4 illustrates the geometric space for this study as sketches of cross sections through the 
crack plane arranged in terms of crack depth-to-thickness ratio (a/B) and crack depth-to-half-length 
ratio {ale) with 0.2 <a/c< 1 and 0.2<al B<().H for a total of 20 different geometries. For each a!B and 
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Figure 4. Geometric space. 


ale combination in figure 4, the smaller, upper illustration is a sketch of the crack plane cross section 
drawn in proportion to the other geometries. (The illustrations for ale- 0.2, a/B = 0.6, and ale = 0.2, 
£/A6 = 0.8 are half-symmetry drawings to allow space for the proportional sketches.) These sketches 
allow the reader to visualize the difference in overall cross section size for each geometry. For each 
a!B and ale combination in figure 4, the lower illustration is a close-up view of the crack plane cross 
section with the thickness held constant for all geometries. The close-up sketches allow one to bet- 
ter see the semielliptical crack shape in relation to the specimen thickness. For all of the geometries, 
B- 1 and L/W= 2. Figure 4 lists the 2c, W, and L values for all the geometries. The plate widths were 
set equal to the greater of W- 5 * 2c or W- 5 * B to minimize width effects on the /-integral solutions 
and to ensure that the plates maintained a ‘plate like’ width-to-thickness aspect ratio for small cracks. 
Utilizing these minimum width criteria precludes the need to include the W/2c ratio as a third vari- 
able in the geometric space. The rationale for the width criteria is discussed in appendix A. 
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2.3 Finite Element Models 


A total of 600 nonlinear finite element analyses were required to perform the analysis of the 
30 material and 20 geometric combinations. All of the finite element models were created using the 
commercial finite element mesh creation and post-processing tool, FEACrack, 1 and the analyses 
were performed using the freely available research code, WARP3D version 16. 3.1. 18 

2.3.1 Finite Element Model Details 

All of the surface cracked plates are modeled with 3D quarter-symmetric FEMs using 
20-node reduced integration isoparametric elements (element type q3disop in WARP3D). A typi- 
cal FEM mesh is shown in figure 5. The number of elements range from 4,520 to 26,928, and the 
number of nodes vary from 21,332 to 122,396 as listed in figure 4. The crack fronts are surrounded 
with ten rings of elements forming the domains for calculation of the /-integral. The innermost 
ring of elements incident on the crack front contains 20-node hexahedrons collapsed into wedges 
with initially coincident nodes along the crack front left unconstrained to permit blunting deforma- 
tions. The number of nodes in the (j) direction along the crack front varies from 49 to 217 as listed in 
figure 4 (CF nodes). The /-integral is evaluated numerically using the domain integral method as 
implemented in WARP3D over each of the ten concentric element domains at each crack front node 
location. The result at the outermost domain (domain 10) is used for the /-integral values in this 
study, and the domain convergence is discussed later in section 3. The reported CMOD values are the 
displacements of the node on the crack mouth on the outer surface centerline of the plate as shown in 
figure 5. 
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2.3.2 Boundary Conditions and Loading 

The FEMs are loaded with 20 to 30 uniform load steps with an average of 2 to 5 Newton 
iterations for convergence within each step to a tight tolerance on residual nodal forces. Symmetry 
planes are enforced through setting the normal nodal displacements on the plane equal to zero. Uni- 
form axial displacements, <5 far , are applied to all of the nodes on the top surface of the plate to apply 
tension. The intent was to apply sufficient displacement to deform the models far into an elastic- 
plastic regime, but not so much that the models would fail to converge on the last analysis step. Since 
hundreds of nonlinear models were run, it was not practical to ‘hand pick’ a unique value for each 
displacement boundary condition — an automated method had to be employed. 

The b far boundary conditions for the FEMs are applied in a two-step process. First, the 
LPPL equation is used to calculate the net section stress, cr net , corresponding to e /;/ = 0.1%. This 
stress is then used to calculate a linear approximation of the far-held displacement using the 
equation 


$ - °net L 

ffiar “ 


2 E 


( 2 ) 


where the factor of Vi is due to symmetry. All of the FEMs were successfully run with the <5far values 
from equation (2) applied, and the models were checked for the deformation level, M, where M is 
defined as 


M = 


j A 


(3) 


In equation (3) .J 0 is the /-integral value at a given parametric crack front angle, b, and is 
the characteristic length for a given 0. For the surface crack problem, two characteristic lengths at 
any 0 location are defined as r^ a and r^ b as shown in figure 6. Considering a 2D slice through the 
specimen along a path defined by r^ a and r^, then r^ a is analogous to a crack length measurement 
and r^ b relates to a remaining ligament measurement. 



Figure 6. Cross section through the crack plane illustrating the 
characteristic lengths, r ( ^ a and r^ h . 
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The authors desired that M~ 20 or less at the final load step (note that M decreases as defor- 
mation increases). The M- values were checked at </> = 30° and d = 90° for both and r^ h for all 
FEMs. For the next step in refining the 5 far boundary condition, the FEMs that did not reach the 
desired deformation level are identified, and their M and <5 tar data are used to predict a new 5 far 
required to reach M~ 20. The method for predicting the new 5 far is illustrated for one of the FEMs 
(a/B = 0.4, ale = 0.6, n = 6, and Elo = 500) in figure 7 wherein log 10 (M) is plotted against the normal- 
ized far-held displacements, 25 far IL, at 0= 30° and 0= 90° for both r ^ and r^ h . In figure 7, deforma- 
tion increases from right to left in the plot. A straight line is fit to the last three data points in the 
result set, and the <5far required to reach M= 20 is estimated for each location from the linear predic- 
tion. The minimum estimated 5 far for each model was applied as a new boundary condition, and this 
subset of FEMs was reanalyzed. The final set of solutions was then collected into a Matlab 19 data 
structure for use in the interpolation routines. 



Figure 7. Method to choose far-held displacement values for FEMs based 
on a desired deformation level, M. 


2.4 Automation Methods 

Handling this large array of nonlinear models was made practical by computer routines that 
automate the process of building the FEMs, running the nonlinear analyses, post-processing model 
results, and organizing the solution results into multidimensional arrays. Computer routines were 
written in Matlab to create the hie storage directory structure and serve as the overall controller for 
the model building, execution, and post-processing procedure. 
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The FEMs were built using FEACrack in batch control mode on a 64-bit Windows XP com- 
puter. FEACrack uses a text-based input hie to control all the model parameters for a given FEM. 
The authors developed Matlab scripts that read in a ‘baseline’ FEACrack input hie for a surface 
crack plate and modified the hie to include the desired plate dimensions, crack dimensions, mate- 
rial properties, and boundary conditions. The script then executes FEACrack in batch mode to 
build the appropriate input hie for WARP3D and saves the hie with an appropriate filename in the 
desired directory. This process was automated to prepare the required 600 input hies for analysis. The 
WARP3D input hies were then moved to a 64-bit Linux server for efficient parallel processing analy- 
sis. Once the FEAs were complete, the WARP3D binary packet result hies were transferred to the 
64-bit Windows XP computer and were post-processed in batch mode with FEACrack using another 
set of Matlab generated scripts. These post-processing steps resulted in a set of 600 text-based result 
hies from FEACrack that contained all of the pertinent model result data. A set of Matlab scripts 
were then used to consolidate the full data set into arrays of /-integral versus 0 values, far-held 
stresses, and CMOD values in an easily indexed Matlab data structure. 

2.5 Interpolation Methodology 

The interpolation methods developed and implemented in this study involve a normalization 
scheme used to scale the solutions. This scheme is described below along with a brief exemplar case 
study. The interpolation process takes place in this dimensionless, normalized space. The interpola- 
tion methods are described on a step-by-step basis. 

2.5.1 Normalization Scheme 


To derive useful results from the solution space, interpolation within the geometry and mate- 
rial dimensions is necessary, but scaling of the solutions with respect to geometry and material is also 
required. To simplify scaling of the solution space, it is advantageous, but not required, to have the 
solution space normalized to a dimensionless state. In this fashion, once interpolation is complete 
within the space, the resulting dimensionless interpolated solution can be easily scaled by to a dimen- 
sioned state by multiplying by representative length and stress scaling factors for the actual geometry. 
There are three primary results in the solution set that need to be normalized: /, CMOD , and far- 
held stress, a. By dimensional analysis, it is clear that the /-integral is conveniently normalized by 
a product of stress and length, therefore the normalized /-integral value, J n , can be written as 


/ 


n 


j 


°ys B 


(4) 


In this case, the normalizing stress is chosen to be the material yield strength (proportional 
limit) represented in the LPPL stress-strain curve model, and the normalizing length is chosen to be 
the plate thickness, B. These are particularly convenient normalizing factors because, as discussed 
previously, both o and B were defined to have unit value in the model space. Thus, the /-integral 
value from the analysis does not change when normalized. The same follows for the CMOD and 
far-held stress results where 


CMOD n 


CMOD 


(5) 
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(7 


(6) 


and the normalized far-held stress, o n , is 


(7 


n 



Starting from this normalized space, a dimensional result is easily obtained by reversing the 
normalization by multiplying by the appropriate dimensional factors, B D and o D , the thickness 
and yield strength, respectively, of the actual specimen or structure. Letting J D , CMOD d , and o D 
represent the dimensional results; 


J D =Jr 


' ysD 


B 


D 


(7) 


CMOD D = CMOD n * B d , 

and 


°D = °n 


' ysD 


( 8 ) 

(9) 


The dimensional force, P D> associated with the far stress, a D , can be obtained by simple 
mechanics following similar form: 


l D 


w d *b 


D 


( 10 ) 


where W D is the dimensional width of the plate. 

2.5.2 Normalization Case Study 

To illustrate the ability to scale nonlinear surface crack solutions across proportional 
geometry changes and differing yield strength values, two surface crack FEMs were analyzed. The 
hrst, model 1, resides near the center of the geometric and material dimensions of the solution space: 
a!B = 0.6, a/c= 0.6, n = 6, E/a ys = 300 with v=0.3. The normalizing dimensional factors for model 1 
are also set to render a ‘normalized’ solution: a, = 1, B= 1. As was the case for other models in the 
space, W is set to 5x2 c, so W- 10, and L = 2x W, so L- 20. 

The second model, model 2, has a proportional geometric scale three times model 1 (B= 3, 
2c = 6, W= 30, L = 60). The material properties for model 2 ( Eto vs , v, n ) are also identical to model 1 
with the exception of <J ys = 50. Note that since Elo ys = 300, the modulus in model 2 is also scaled 
with o vs . Any choice of consistent units may be assumed, and the authors chose the SI units 
(mm for length and MPa for stress) for this example. Figure 8 shows a half-symmetric sketch of 
the two models drawn in proportion to one another. Model 2 was run in displacement control to 
three times the far-held displacement, 5 far of model 1 . The model 1 results were then scaled using 
equations (7)— (10) with o ysD = 50, B D = 3, and 1^=30 to predict the model 2 result. Model 2 and 
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Model 1 


W 



Figure 8. Cross section through the crack plane for two surface crack geometries 
with different physical dimensions but the same geometric proportions. 


scaled model 1 results for P versus CMOD, / total , and J e/ (linear-elastic portion of J total ) versus 
CMOD, and / total and J e j versus (f) at the last load step are compared in figures 9-1 1 illustrating per- 
fect agreement for each variable. 



Figure 9. Comparison of the force versus CMOD response for models 1 and 2. 
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Figure 10. Comparison of / total and J e j versus CMOD for models 1 and 2. 



Figure 11. Comparison of / total and J e j versus <p for models 1 and 2. 
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2.6 Solution Space Interpolation 


Having established the basis for the nondimensional solution space and procedures for scal- 
ing proportional geometry and material yield strength to dimensioned values, consideration turns to 
interpolation within the space, which provides an estimated solution at any crack shape and depth 
within the geometric space and at any modulus of elasticity and strain hardening exponent within 
the material space. The quantity of interest is the J n value as a function of 0 around the crack perim- 
eter, J n {<p). For each of the 600 models in the space, J n (<p) is calculated as a function of increasing 
load. The state of the load can be described by either the models’ far-held stress, o n , or displace- 
ment at the crack mouth, CMOD n . Though o n (or force) is an intuitive descriptor of the load, for 
elastic-plastic analysis, the CMOD is a more reliable predictor of J ( J is nearly a linear function of 
CMOD in the plastic regime); 7 ’ 17 ’ 20 so the authors chose to use CMOD as the quantifier of the load- 
ing condition in the interpolation methodology. Figure 12 shows a plot of / n (0) versus CMOD n for 
a solution with 30 load increments. Open symbols are placed at the 0 = 30° location to help visualize 
the J n versus CMOD trajectory for a given (j) location. In the solution space, the relationship between 
a n and CMOD n is maintained, thus by dividing a n by the final dimensioned area, a prediction of the 
P versus CMOD trace is available. This trace is particularly useful for interpretation of experimental 
surface crack test results. 



Figure 12. J((j>) versus CMOD u space. 
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The solution space consists of 600 result data sets, each containing J n ((j)) versus CMOD n 
and o n versus CMOD u data. The space is structured in a four-dimensional array that is most easily 
visualized by considering a 4x5 geometry matrix with four rows of alB ratios and five columns of 
ale ratios. Within each of the 20 geometric combinations, there exist 30 material solutions described 
by a 5 x 6 matrix of material solutions, five values of n and six values of E/a ys . The solution space 
is readily indexed by these four dimensions. For a given model result, R, the solution is given by the 
notation: R(a/B , ale, n, Elo ). Figure 13 shows a conceptual illustration of the R(a/B , ale, n, Elo ys ) 
solution space with the geometric space at the highest level and the entire material space existing at 
the next level repeated within each geometric combination followed by the J n ((j)) versus CMOD n and 
a n versus CMOD n data for each of the 600 models at the lowest level. 

In general, the actual surface crack geometry and material of interest will not fall 
directly on an existing solution and interpolation is necessary. To interpolate to a new solution, 
R[al B, a I c, n, El <7^J,the first step is to identify the subset of the 600 model space that will be 
active in the interpolation process by determining the location of R in the geometry and material 
matrices. For illustration, consider a choice of / B = 0.5 , a lc = 0.5 , n = 8 , E / <y ys = 400 j that is 
located between the cells labeled gj-g 4 in the geometry matrix of figure 13. The four ‘nearest neigh- 
bor’ subset solutions are the geometry combinations designated as g\-g^. For each of the g, -g 4 
geometries, a point for R can be placed in the material matrix resulting in materials m 1 -m 4 . Identify- 
ing the sets g\-g^ and their associated m |-m 4 sets provides the 16 nearest neighbor data sets for use 
in the interpolation of the R solution. 

The 16 subsets represent different geometry and material combinations, so there is no expec- 
tation that the data sets will have the same nodal 0 spacing, number of load increments, spacing of 
load increments, or final deformation values. Therefore, to simplify the interpolation process, the 0 
locations and load increments are adjusted to be uniform across the 16 subset solutions. Even spac- 
ing in 0 is easily created by linear interpolation of the existing J n ((p) data at standard (J) locations for 
each load increment in every model. In this study, the J n values are calculated in 2° increments for 
0° < 0< 90° for a total of 46 J n {^) values at each load increment. 
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Crack Depth-to-Thickness Ratio, a/6 



Figure 13. Conceptual illustration of the interpolation space. 
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The methodology for setting uniform spacing and magnitude for loading conditions across 
the 16 subset solutions is not as straightforward. Since the FEMs are run to different levels of o n 
and CMOD n , several choices must be made concerning how to divide the solutions into load incre- 
ments including: the measure used for parsing (a n or CMOD n ), the load increment spacing, and 
a method to either truncate or extend the solutions. As stated above, the authors chose to use CMOD n 
as the parsing measure of loading condition, thus a common maximum value of CMOD n across all 
16 sets is needed. Figure 14(a) plots a n versus CMOD n for the 16 subset solutions and illustrates the 
challenge in determining the maximum allowed CMOD n value for the solution set. The minimum, 
maximum, and average CMOD n values are designated by vertical lines in figure 14(a). Truncating all 
of the data sets to the minimum CMOD n of the 16 sets severely limits the extent of the solutions. In 
turn, extrapolating all of the data sets out to the maximum CMOD n value requires extending some 
solutions to CMOD n values far beyond their final converged load increment. Instead, the authors 
chose to set the maximum common loading condition equal to the average CMOD n of the 16 data 
sets. Solutions that extend beyond the average CMOD n are truncated. Solutions with maximum 
CMOD n values less than the average are extrapolated out to the average CMOD n by fitting a power 
law of the form o n =a* (CMOD n )P+yto the last five o n versus CMOD n data points in the set, where 
a, [3, and yare fitting constants. Once the 16 data sets have all been either truncated or extrapolated 
to the same maximum CMOD n , each set is divided into 20 even CMOD n increments common across 
all sets. Linear interpolation is used to determine the values of a n for the 20 CMOD n increments. 
Figure 14(b) illustrates the 16 a n versus CMOD n data sets with 20 even CMOD n increments out to 
the average CMOD n . 
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Figure 14. Illustration of (a) effect of truncating the solution space based 
on the average of the maximum CMOD n values on the a n versus 
CMOD n results, and (b) dividing the CMOD n into even incre- 
ments for a n versus CMOD n interpolation. 
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With each of the 16 subset solutions adjusted to common increments of CMOD n , the J n ( (p) 
values for the new load increments must be interpolated. Figure 15(a) shows J n ((j)= 90°) versus 
CMOD n for the same 16 subsets shown in figure 14. The J n versus CMOD n solutions that extended 
beyond the average CMOD n are truncated while the solutions that are less than the average CMOD n 
are extrapolated using linear extrapolation. Linear extrapolation works well for these data because 
the J n versus CMOD n response becomes approximately linear once into the elastic-plastic regime. 
Figure 15(b) illustrates J n ((t> = 90°) versus CMOD n data sets with 20 even CMOD n increments out to 
the average CMOD n . 


21 



Figure 15. Illustration of (a) effect of truncating the solution space based 
on the average of the maximum CMOD n values on the J n versus 
CMOD n results, and (b) dividing the CMOD n into even incre- 
ments for J n versus CMOD n interpolation. 
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Once the 16 subset solutions are prepared with even 0 spacing and load increments, the 
step-by-step interpolation process to the final R can proceed. Figure 16 illustrates the four-step 
interpolation process used by the authors. Each of the steps are indicated by the Roman numerals 
I through IV. At each interpolation step I through IV, J n is interpolated at each 0 location for each 
CMOD n increment to determine a new .J n ( d) versus CMOD n . A second interpolation process evalu- 
ates o n for each CMOD n increment. Step I in the interpolation process begins with the 16 prepared 
subset solutions for linear interpolation based on crack shape ratio, ale. Within the 16 solutions, 
there are eight pairs of solutions that differ only by ale. Linear interpolation occurs between these 
pairs based on the target ale value, resulting in eight new solutions, all with the same ale value but 
different combinations of a/B, n, and Elo ys . In step II, linear interpolation is performed on the four 
pairs of solutions that differ only by crack depth, alB , to the chosen a!B value for R . This results 
in four new solutions with the same geometric values but different material values. In step III, two 
pairs of solutions are linearly interpolated in terms of n to the chosen n value for R resulting in two 
new solutions differing only in their E/o vs values. Finally, in step IV, the remaining two solutions are 
interpolated to the chosen E/a ys , resulting in the final R interpolated solution. As will be shown in 
section 3, the J „(<!>) and a n values have nonlinear spacing with Elo ys ; therefore, for the step IV inter- 
polation, linear interpolation was performed in terms of log 10 (F/(J vs> ) and the log 10 of the chosen 
E/a,,„ for R value. 

A * 3 
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Figure 16. Conceptual illustration of the interpolation process. 
Arrow indicates interpolation direction. 


Figures 17-20 illustrate the interpolation steps I through IV, respectively, for the example 
interpolated solution, / B,a / c,n,E / o ys j = ^(0.5,0.5,8,400) . Each figure shows a comparison of 
results for: (1) o n versus CMOD n , (2) J n (<j> = 90°) versus CMOD n , and (3) for the final common 
loading condition of the nearest neighbor solution subsets along with the final interpolated result, 
R . The solution subset results are plotted as solid lines with open symbols; the open symbols repre- 
sent specific materials (m|-m 4 in figs. 13 and 16) while the colors represent specific geometry choices 
( gi~g 4 in figs. 13 and 16). The final interpolated results are plotted with a solid line with filled star- 
shaped symbols. Figure 17 represents the original 16 subsets prior to step I interpolation process on 
ale. Figure 18 shows the resulting eight solutions of step I prior to the step II interpolation for alB. 
Figure 19 shows the resulting four solutions of step II prior to step III interpolation on the n value. 
Lastly, figure 20 shows the final two solutions to be interpolated by Elc J in step IV. 
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Interp., alB= 0.5, alc= 0.5, n= 8, El o =400 


. a/6=0.4, a/c=0.4, n=6, £/cr ys = 300 ~| 
■ a/8=0.4, a/c= 0.6, a=6, £/cr ys =300 
. alB= 0.6, a/c=0.4, n=6, Elo y =300 
. a/6= 0.6, a/c= 0.6, n=6, £/cr ys =300 __ 


■ alB= 0.4, a/c=0.4, n=6, Elo ys =500 

■ a/8=0.4, a/c= 0.6, n=6, Elo y =500 

■ alB= 0.6, a/c=0.4, n=6, Elc ys =500 

■ alB= 0.6, a/c= 0.6, n=6, Elc ys =500 

. alB=0A, a/c=0.4, n=10, £/cr ys =300 

■ a/B=0.4, a/c= 0.6, n=10, Elo y =300 
. a/B=0.6, a/c=0.4, n=10, Elo y =300 

■ alB= 0.6, a/c= 0.6, n=1 0, Elo y =300 
. a/B=0.4, a/c=0.4, n=10, £/cr ys =500 
. a/B=0.4, a/c=0.6, n=10, Elo ys =500 
. a/B=0.6, a/c=0.4, n=10, £/cr ys =500 

■ a/B=0.6, a/c=0.6, n=10, £/cr ys =500’ 




(b) 


CMOD n 


(c) 


Parametric Angle, 0 (deg) 


Figure 17. Comparison of results for (a) cr ;7 versus CMOD n , (b) /,/</) = 90°) versus CMOD u , 
and (c) J n versus d for the final common loading condition for the 16 solution 
subsets and the final interpolated result prior to interpolation step I to the desired 
ale value. 
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- Interp., a/B= 0.5, a/c= 0.5, n= 8, Ela ys = 400 


- a/8=0.4, a/c= 0.5, n= 6, £/cy s = 300 

— e— 

- a/B= 0.6, a/c= 0.5, n=6, £/cr ys =300 

B— 

a/B=0.4, a/c= 0.5, n=6, Ela ys =500 

B— 

alB= 0.6, a/c= 0.5, n=6, £/cr ys =500 

— e- 

- a/B=0.4, a/c=0.5, n=10, £/cr ys =300 

— e- 

alB= 0.6, a/c= 0.5, n=10, £/cr ys = 300 

-A- 

- a/B=0.4, a/c=0.5, n=10, Elo y =500 

—A- 

- alB= 0.6, a/c= 0.5, n=10, Elc y s =500 



(b) CMOD n 



(c) Parametric Angle, <f> (deg) 


Figure 18. Comparison of results for (a) o n versus CMOD n , (b) J n ((j) = 90°) versus CMOD n , 
and (c) J n versus (j) for the final common loading condition for the eight solu- 
tion subsets and the final interpolated result prior to interpolation step II to the 
desired a!B value. 
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- Interp., alB= 0.5, a/c= 0.5, n= 8, £/cr ys =400 


- alB= 0.5, alc= 0.5, n=6, Elo ys = 300 

B— 

alB= 0.5, alc= 0.5, n= 6, Ela y =500 

— 0- 

- alB= 0.5, alc= 0.5, n=10, Elo y =300 

—A- 

alB= 0.5, a/c=0.5, n=10, Ela ys =500 



(b) CMOD n 



Parametric Angle, (f> (deg) 


Figure 19. Comparison of results for (a) o n versus CMOD n , (b) .J n ((j) = 90°) versus 
CMOD n , and (c) J n versus 0 for the final common loading condition 
for the four solution subsets and the final interpolated result prior to 
interpolation step III to the desired n value. 
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Interp., alB= 0.5, alc= 0.5, n= 8, £/cr ys =400 


-e — alB= 0.5, alc= 0.5, n=8, £/cr ys =300 
-a — alB= 0.5, alc= 0.5, n=8, £/cj ys =500 




Figure 20. Comparison of results for (a) o n versus CMOD n , (b) .J n ((j) = 90°) versus 
CMOD n , and (c) J n versus 0 for the final common loading condition 
for the two solution subsets and the final interpolated result prior to 
interpolation step IV to the desired Elo s value. 


29 


3. REVIEW AND VERIFICATION 


The solution space described herein is too broad for detailed review; however, a brief para- 
metric review of the space is warranted to convey to the reader the influence of the four parameters 
of the space. Following this brief review, the authors present the results of the verification process 
used to ensure the reliability of the solution space and the interpolation routines. The verification 
process also serves to establish the expected error in the interpolated solutions relative to direct 
modeling results. The verification process included four steps: (1) a review of the linear-elastic stress 
intensity solutions at the first model load increment, (2) confirmation of elastic-plastic model con- 
vergence through domain integral independence at the final model load increment, (3) comparison 
of interpolated solutions against a 25 model benchmark solution set, and (4) review of an existing 
round robin (RR) on elastic-plastic surface crack analysis and a comparison of an interpolated solu- 
tion relative to the population of independent RR solutions. 

3.1 Representative Model Results 

The following parametric study provides a brief review to help visualize the influence of the 
parameters in the interpolation space. The same model used in section 2.5.2 with a/B = 0.6, alc = 0.6, 
n = 6, and Elo ys = 300 was chosen as a representative result set from near the center of the geometry 
and material matrices. As discussed in the case study, these parameters fall directly on a model within 
the space requiring no interpolation. Each of the parametric variations described also fall directly on 
existing solutions. By holding the representative model parameters constant and individually varying 
a/B , ale, n, and Elo ys across their respective matrix values, useful figures can be made demonstrating 
the effects of varying geometry and material values. Each of the figures 21-24 has three subplots 
for the representative model and various choices of alB, ale, n, and E/o : (a) o n versus CMOD n , 
(b) ./„(</) = 90°) versus CMOD n , and (c) J n ((j)). In figures 21-24, the maximum common CMOD n value 
for all of the FEMs being compared is designated with a vertical line in parts (a) and (b) of the fig- 
ures. The J n ((t>) results plotted in part (c) correspond to the maximum common CMOD n . It is also 
important to note that the results presented in figures 21-24 are data directly from the FEM; they 
have not yet been interpolated into uniform (j) or load increments. 
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Figure 21 illustrates the effect of varying crack depth, alB, with ale, n, and F/c^.held constant 
for the representative geometry. Changes in alB have a nearly uniform effect on a n versus CMOD n , 
J n ((p= 90°) versus CMOD n , and ./,/</)). In figure 21(c), the magnitude and shape of the J n -(j) curves 
differ primarily because at the maximum common CMOD n point, the a!B = 0.2 and 0.4 analyses are 
nonlinear in their global response, while the other a!B results are still linear-elastic. 



— e— 

-a/e = 0.2 

B— 

-a/e = 0.4 

— 0- 

- a/6 = 0.6 

— A — 

-a/6 = 0.8 



Figure 21. Comparison of results for (a) o n versus CMOD n , (b) J n ( (j) = 90°) versus 
CMOD n , and (c) versus (h for various aJB values, ale = 0.6, n = 6, and 
F/(j = 300. 

S'* 
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Figure 22 illustrates the effect of varying crack aspect ratio, ale, with alB , n, and EIg held 
constant. Changing ale does not have a uniform effect on o n versus CMOD n . There is a significantly 
different response for the ate = 0.2 and ale- 0.4, a n versus CMOD n curves indicating that perhaps 
a set of results at ale = 0.3 would be useful to improve interpolation between these curves. Parts (b) 
and (c) of figure 22 confirm that J n versus CMOD u at (j)= 90° is controlled primarily by crack depth, 
a/B. The crack aspect ratio drives J n variations for other (j) locations. 



— e— 

- a/c= 0.2 

B— 

- a/c= 0.4 

— 0- 

a/c= 0.6 

— A — 

a/c- 0.8 


a/c= 1 


0 0.02 0.04 0.06 



Figure 22. Comparison of results for (a) a n versus CMOD n , (b) .J n ((j) = 90°) versus CMOD n , 
and (c) J n versus (j) for a/B =0.6, various ale values, n = 6, and E/o vs = 300. 


32 


Figure 23 illustrates the effect of varying strain hardening exponent, n, with alB. ale, and 
E/o ys held constant. As expected, changes in n have no effect on the linear portion of the o n versus 
CMOD n response but do affect the curve once nonlinear. Following directly from the LPPL tensile 
response, the n = 3 model has the highest o n for a given CMOD n while the n = 20 nonlinear response is 
almost flat and has the lowest o n for a given CMOD n . In figure 23(b), the J n ((j) = 90°) versus CMOD n 
results are identical for all n values for low deformations in the LEFM regime and then begin to 
diverge for higher deformations as plasticity increases. Compared to the other parameters, the strain 
hardening exponent has the weakest effect on the solution space. 



Figure 23. Comparison of results for (a) a n versus CMOD n , (b) ./„(</) = 90°) versus CMOD n , 
and (c) J n versus b for a!B = 0.6, various n values, and E/o = 300. 
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Figure 24 illustrates the effect of varying the reciprocal of the yield strain, E/<j ys = 1 le with 
a/B, ate, and n held constant. In figure 24(a), the a n versus CMOD n response is very sensitive to 
the Elo ys value and is not a linear function of Elo ys . The ./ ;) -d results in figure 24(c) also appear to 
be a nonlinear function of Elo ys . These observations led to the choice of linear interpolation for 
Elo being performed in terms of log 10 (£yc7 ra ) as discussed previously in section 2.5. Interpolation 
in terms of e would also be a potential option for future investigation. These observations also 
suggest that a set of results at E/a vs = 150 would be useful to improve interpolation between the 
E!o vs = 1 00 and 200 solutions. 
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Figure 24. Comparison of results for (a) a n versus CMOD n , (b) J n ((j) = 90°) versus CMOD n , 
and (c) J n versus 0 for a/B =0.6, ale = 0.6, n = 6, and various E/o ys values. 


3.2 Verification Process 


3.2.1 Linear-Elastic Result Comparisons 

The linear-elastic J n results were compared with linear-elastic stress intensity results 
calculated using Newman-Raju equations. 4 This comparison confirms the basic quality of 
the 20 FEM geometries used in this study by reaffirming the well-established NR solutions for 
surface cracks in tension. To establish a basis for comparison, the J n ((\j) results from the first load 
increment for each of the 20 geometries were converted to the mode I linear-elastic stress intensity 
factor, Kj, through the plane strain equation: 



( 11 ) 


The specific material combination used for this comparison is not important because Kj is 
independent of material properties, and J n can be considered completely elastic at the first load 
increment in each analysis. A geometry correction factor, A, was then calculated for each of the 
20 configurations using the equation 


A = 



( 12 ) 


The specific stress and geometry of each FEM was input into the NR equations to calculate 
A NR for each 0 location along the crack front using equation (12). Figure 25 shows the results of 
the comparison between the linear-elastic FEM results and the NR equations by plotting the ratio 
Afem^NR versus <f>. A solid, heavy horizontal line is drawn at a ratio value of 1, indicating perfect 
agreement for all (f) locations, and dashed horizontal lines are drawn above and below to indicate 5% 
deviation. The symbols in figure 25 represent alB values and the colors represent ale values. Newman 
and Raju 3 report that their surface crack stress intensity equation was within ±5% of their FEM 
results for a/B< 0.8, but the accuracy of the equation is unknown for alB> 0.8. Overall, the FEMs 
and the NR solutions are in excellent agreement with only two models with a/B=0.S deviating by 
more than 5% at the depth location for deep cracks. In general the A FEM values are marginally less 
than the A NR values, confirming the slightly conservative nature of the NR equations for Kj calcula- 
tions. 
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Figure 25. Comparison of linear-elastic stress intensity correction factors from the 

NR equations to the values calculated from the 20 geometries in this study. 


3.2.2 J-Integral Solution Convergence Checks 

The next step in the verification process is to confirm the quality of the elastic-plastic /-inte- 
gral solution at the final load increment of each model as indicated by the path independence of 
the /-integral. All 600 FEMs used in the study are confirmed to have converged /-integral values 
based on this metric, /-integral convergence becomes more difficult with increasing plasticity, thus 
checking convergence at the final load increment confirms the quality of the complete analysis. 
Figure 26 demonstrates the path dependence of the J n values at 0 = 90° for the last analysis step in all 
600 models. In this figure, the horizontal axis is the domain integral number with 1 being the inner- 
most domain nearest the crack tip and 10 being the outermost domain. The vertical axis represents the 
J n values at 0= 90° for the last load increment normalized by the maximum J n domain value for that 
FEM at 0=90°; therefore, each domain location in figure 26 contains 600 data points. As expected, 
the largest deviation from the outermost domain value is nearest the crack tip at domain 1, and the J n 
values rapidly converge towards the outermost domain value. As a check of domain convergence, the 
percent difference between the domain 9 and 10 J n values is calculated for every 0 location in all the 
models at the last load increment. Figure 27 illustrates the percent difference between domain 9 and 10 
J n values at every (f) across the entire solution space. The FEMs demonstrate excellent J n convergence 
with most of the data having 0.5% or less difference between the last two domains. Slightly higher dif- 
ferences exist at the free surface location, but the differences still remain below 0.85%. From this plot, 
one concludes that the FEMs all had a sufficient number of domains and a large enough outer domain 
size to ensure sufficiently converged J n values at the last load increment of every analysis. 
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Figure 26. Normalized J n (<j)= 90°) for each domain at the final load step 
for all 600 solutions. 
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Figure 27. Percent difference in domains 9 and 10 J n values for all 
(f) locations at the final load step for all 600 solutions. 
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3.3 Elastic-Plastic Solution Comparison to Benchmark Data Set 
3.3.1 Description of Benchmark Data Set 

The quality of the linear-elastic FEM solutions was confirmed through comparison to the 
Newman-Raju equations. Unfortunately, no such established solution exists for comparison of 
the elastic-plastic /-integral results; therefore, the authors developed a set of elastic-plastic ‘bench- 
mark’ solutions consisting of 25 FEMs to verify the surface crack solution interpolation method- 
ologies. Five surface crack geometries and five materials were chosen to comprise the benchmark 
solutions set. The geometry and material choices are given in table 1 and were chosen to fall in the 
‘gaps’ in the geometry and material matrices. The benchmark FEMs were built, executed, and post- 
processed using the same methods described in section 2. All of the benchmark FEM results were 
then collected in a Matlab data structure. 


Table 1 . Benchmark FEM geometry and material values. 


Geometry Values 

a/6 

ale 

0.3 

0.3 

0.3 

0.9 

0.5 

0.5 

0.7 

0.3 

0.7 

0.9 


Material Values 

n 

ElCys 

3.5 

150 

3.5 

850 

8 

400 

15 

150 

15 

850 


3.3.2 Comparison of Interpolated Solutions and Benchmark Dataset 

A Matlab script was written to evaluate an interpolated solution for each of the 25 bench- 
mark solutions. The script then creates a series of plots and tables comparing the benchmark and 
interpolated solutions to aid in investigating the effectiveness of the interpolation methodology. 
Table 2 lists the ale , alB , n, and Elo s values for the benchmark models and assigns a ‘model No.’ 
label to each. For each benchmark case, two values of CMOD n ( CMOD n A and B) were cho- 
sen as displacement levels at which to compare the results. CMOD n A represents a point slightly 
past the LEFM regime, just beyond the rollover in the P versus CMOD n trace. CMOD n B is the 
point of highest common load increment between the benchmark FEM and the interpolated solu- 
tion and represents the EPFM regime with large plasticity. The regime between CMOD n A and 
CMOD n B is the intended range for use of this solution set. Table 2 contains the percent difference 
between the benchmark solutions and the interpolated solutions for P, J n ((t>= 30°), and J n ((j> = 90°) at 
CMOD n A and B. Overall, the interpolated solutions are a nice match to the benchmark solutions. At 
CMOD n A, the average percent difference for all compared values is less than 5%, and at CMOD n B 
the average percent difference is less than 3%. The largest percent difference at CMOD n A is 1 1% for 
the model No. 5 J „((/>= 30°) comparison, and the largest percent difference at CMOD n B is 6% for the 
model No. 19 J n ((j> = 30°) comparison. 
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Table 2. Benchmark model numbers, model descriptions, and interpolated result comparisons. 



% Difference at CMOD A 

% Difference at CMOD B 

Model No. 

alB 

ale 

n 

Ho vs 

P 

J n (0 = 3O°) 

J n U>= 90°) 

P 

J n (0 = 3O°) 

J„(0=9O°) 

1 

0.3 

0.3 

3.5 

150 

4.4 

5.5 

1.8 

- 2.4 

4.5 

3.9 

2 

0.3 

0.3 

3.5 

850 

- 1.0 

3.6 

0.8 

- 3.7 

1.5 

1.5 

3 

0.3 

0.3 

8.0 

400 

7.5 

9.2 

3.9 

- 1.9 

3.2 

3.4 

4 

0.3 

0.3 

15.0 

150 

9.3 

8.2 

4.2 

2.1 

3.7 

4.9 

5 

0.3 

0.3 

15.0 

850 

6.8 

11.0 

3.9 

- 1.0 

4.1 

4.5 

6 

0.3 

0.9 

3.5 

150 

1.7 

1.9 

0.7 

- 3.3 

1.4 

1.3 

7 

0.3 

0.9 

3.5 

850 

- 0.4 

1.4 

0.5 

- 3.6 

- 0.2 

0.1 

8 

0.3 

0.9 

8.0 

400 

1.6 

5.1 

3.9 

- 1.7 

1.5 

1.5 

9 

0.3 

0.9 

15.0 

150 

7.5 

5.8 

4.2 

0.2 

4.0 

3.7 

10 

0.3 

0.9 

15.0 

850 

3.8 

6.5 

5.4 

- 1.0 

2.3 

2.0 

11 

0.5 

0.5 

3.5 

150 

2.1 

3.5 

1.2 

- 1.8 

2.8 

2.0 

12 

0.5 

0.5 

3.5 

850 

- 1.0 

1.8 

0.0 

- 2.0 

0.5 

0.4 

13 

0.5 

0.5 

8.0 

400 

3.1 

4.8 

1.8 

- 1.7 

2.0 

1.0 

14 

0.5 

0.5 

15.0 

150 

6.7 

4.8 

1.7 

1.5 

4.1 

2.4 

15 

0.5 

0.5 

15.0 

850 

3.4 

4.7 

1.0 

- 1.4 

2.5 

1.1 

16 

0.7 

0.3 

3.5 

150 

4.7 

5.8 

- 0.4 

- 1.6 

5.4 

0.7 

17 

0.7 

0.3 

3.5 

850 

3.0 

4.6 

- 1.3 

- 1.4 

2.5 

- 2.0 

18 

0.7 

0.3 

8.0 

400 

5.2 

6.1 

- 2.3 

0.8 

4.0 

- 1.3 

19 

0.7 

0.3 

15.0 

150 

6.4 

5.9 

- 1.4 

1.4 

6.0 

0.4 

20 

0.7 

0.3 

15.0 

850 

4.2 

5.8 

- 3.0 

0.2 

3.7 

- 1.2 

21 

0.7 

0.9 

3.5 

150 

3.1 

1.9 

0.5 

- 1.9 

0.7 

- 1.1 

22 

0.7 

0.9 

3.5 

850 

1.4 

1.0 

- 0.3 

- 2.0 

- 0.8 

- 2.4 

23 

0.7 

0.9 

8.0 

400 

3.3 

1.7 

- 0.5 

- 1.1 

0.7 

- 1.7 

24 

0.7 

0.9 

15.0 

150 

5.3 

1.9 

0.1 

1.4 

2.9 

1.2 

25 

0.7 

0.9 

15.0 

850 

1.0 

1.4 

- 0.9 

- 0.8 

1.0 

- 0.8 


Max . 

9.3 

11.0 

5.4 

2.1 

6.0 

4.9 

Avg . 

3.7 

4.6 

1.0 

- 1.1 

2.6 

1.0 

Min . 

- 1.0 

1.0 

- 3.0 

- 3.7 

- 0.8 

- 2.4 

Std . Dev . 

2.7 

2.6 

2.2 

1.6 

1.7 

2.0 
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The aforementioned Matlab script creates a series of six comparison plots for each of the 
25 benchmark FEMs. Figures 28 and 29 show two examples of the comparison plot series for model 
No. 1 and model No. 25, and the complete set of comparison plots are included in appendix B. 
Figures 28 and 29 each have six plots (a)-(f) that are described as follows: 



(c) 



Figure 28. Comparison of benchmark model No. 1 FEM results and interpolated results for 
(a) P versus CMOD n , (b) M versus CMOD n , (c) benchmark FEM/interpolated 
results versus CMOD n , (d) J n versus CMOD n at <p= 30° and 90°, (e) J n versus P at 
(p= 30° and 90°, and (f) J n versus <p for CMOD n values A and B. 
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Figure 29. Comparison of benchmark model No. 25 FEM results and interpolated 

results for (a) P versus CMOD n , (b) M versus CMOD n , (c) benchmark FEM/ 
interpolated results versus CMOD n , (d) J n versus CMOD n at (p= 30° and 90°, 

(e) J n versus P at <p= 30° and 90°, and (f) J n versus (p for CMOD n values A 
and B. 

• Plot (a) compares the reaction force, P, versus CMOD n results for the benchmark FEM and cor- 
responding interpolated solution. In this plot, the solid line represents the benchmark FEM solu- 
tion and the red ‘X’ symbols represent the interpolated result at the 20 distinct load increments 
in the interpolated solution. In plots (a)-(d), vertical dashed lines are drawn at CMOD n A and 
CMOD n B. 

• Plot (b) illustrates the deformation level, M, as a function of CMOD n . The four curves represent 

different M - values from equation (3) where r & is r^ a and r^ h at (p= 30° and 90°, respectively. 
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• Plot (c) presents a ratio of benchmark FEM results to interpolated results for P, ./„(</)= 30°), and 
J„((t> = 90°) as a function of CMOD n . Perfect agreement is symbolized by the horizontal line at 
a ratio value of 1 while ±10% variation lines are drawn as horizontal dash-dot lines. 

• Plot (d) illustrates ./„(</) = 30°) and J „((/) = 90°) versus CMOD n for the benchmark FEM and corre- 
sponding interpolated solution. For plots (d)-(f), the solid and dashed lines represent the bench- 
mark FEM solution and the red ‘X’ symbols indicate the interpolated result. 

• Plot (e) shows J n ((t> = 30°) and ./,/</) = 90°) versus P for the benchmark FEM and corresponding 
interpolated solution. Two vertical dashed lines are drawn at the benchmark FEM reaction forces 
corresponding to CMOD n A and CMOD n B. 

• Plot (f) compares J n ( </>) for the benchmark FEM and related interpolated solution for deformation 
levels corresponding to CMOD n A and CMOD n B. 

3.4 Interpolated Solution of the Round Robin Surface Crack Test 

The authors led an interlaboratory round robin concerning the elastic-plastic analysis of sur- 
face cracked plates with 15 participants as documented in NASA/TM — 20 12-2 17456. 17 A surface 
crack tension test was performed to serve as the basis for the RR work. The experiment employed 
a 2219-T8 aluminum specimen with JF=88.82 mm, B- 9.5 mm, L- 177.8 mm (uniform cross sec- 
tion length), a = 6.17 mm, and 2c- 12.7 mm as shown in figure 30. The specimen was loaded under 
displacement control in tension until ductile tearing was detected. The tearing force was 252 kN 
corresponding to a tearing CMOD of 0.1 14 mm. The extent of ductile tearing was then highlighted 
using a combination of sodium hydroxide surface staining and fatigue cycling to mark the surface, 
and the specimen was then loaded to failure. Figure 31 shows a picture of the fracture surface with 
the local location of ductile tearing shown at 0= 17°. 
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Figure 30. Round robin specimen configured for testing. 
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Local Tear Location 
at 0=17° 

Figure 31. Round robin specimen fracture surface with tearing location indicated. 


An interpolated solution was used to analyze the test and this analysis is compared to the 
FEAs of the RR participants. In the following section, the authors’ original FEA performed for the 
RR is labeled ‘FEA,’ and the other participants’ results are labeled ‘labs 2-15.’ An LPPL approxima- 
tion of the material’s stress-strain curve is required to estimate an interpolated solution to the prob- 
lem. The interpolated solutions assume the same elastic properties provided to the RR participants 
(E = 74.46 GPa and v= 0.33) with the exception that v is a fixed value of 0.3 in all of the interpolated 
solutions. Figure 32 illustrates the 2219-T87 engineering stress versus engineering strain test data as 
a dashed line. The RR participants were free to model this material behavior in any fashion. The open 
circle symbols in figure 32 show the authors’ original RR FEA stress-strain table inputs in the form of 
a multilinear approximation of the uniaxial tensile test data. The first data point used in the original 
analysis is at a stress value close to the proportional limit and well below the 0.2% offset yield value. 
Using an LPPL representation of the stress-strain curve does not allow one to place the initial yield 
point this low and still capture the remainder of the plastic behavior, so compromises have to be made 
to attempt a ‘best fit’ of the experimental stress-strain data. For interpolated solutions, the authors 
chose five reasonable fits to the nonlinear portion of the stress-strain curve labeled ‘int. 1’ through 
‘int. 5’ in figure 32, and table 3 lists the a vs and n values assumed for each of the five cases. Int. 1 assumes 
a lower value of a and a lower n value while int. 5 assumes a higher o vs value close to the 0.2% off- 
set yield value and a correspondingly higher n value. All five cases were chosen to represent reason- 
able engineering choices for fitting the material behavior, and the sensitivity to the material choices 
were studied. 


44 


500 



Engineering Strain 


Figure 32. Round robin uniaxial tensile stress versus strain data compared with FEA and LPPL 
material fits. 


Table 3. Round robin FEM and interpolation method sensitivity study result comparisons. 


Result Name 

°ys 

(MPa) 

n 

J(<P= 17°) 
(kJ/m 2 ) 

J % Diff. 
from FEA 

Pat Tearing 
CMOD (kN) 

P % Diff. 
From Test 

P % Diff. 
From FEA 

FEA 

stress-strain table 

19.75 

0.00 

258.9 

2.81 

0.00 

Int. 1 

344.7 

7.5 

19.63 

-0.65 

252.8 

0.40 

-2.35 

lnt.2 

351.6 

8.5 

19.80 

0.22 

255.4 

1.44 

-1.33 

Int. 3 

365.4 

9.5 

20.14 

1.94 

261.6 

3.90 

1.05 

Int. 4 

372.3 

10.5 

20.31 

2.80 

264.3 

4.94 

2.07 

Int. 5 

386.1 

12.5 

20.55 

4.05 

267.6 

6.27 

3.36 
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Figure 33 illustrates the effect of the choice of LPPL material constants on (a) P versus 
CMOD, (b) / versus CMOD at 0= 17°, and (c) ./(</))• The figure compares the interpolated results 
to the authors’ RR FEA results. The results are essentially identical in the linear-elastic regime and 
then diverge as plasticity increases. All of the interpolated solutions extend to a greater deformation 
level than the FEA result. Int. 4 and 5 extend to a higher deformation level than int. 1-3 because 
of the different 16 nearest neighbor subsets used in the solutions due to the differing choices for n. 
Dashed vertical lines are drawn in figures 33(a) and (b) that correspond to the test tearing CMOD 
and the maximum CMOD achieved in the FEA result. Part (c) illustrates the / versus b results at 
these two CMOD levels. The biggest effect of the change in LPPL material properties is in the non- 
linear P-CMOD response; the interpolated /-integral results are fairly insensitive to the given choices 
for plastic material properties. Table 3 gives the /-integral values for </)= 1 7° and the P results at the 
tearing CMOD. The percent difference between the FEA and interpolated solution /-integral results 
are small and range from -0.65% to 4.05%. The percent difference between the FEA and interpolated 
solution force values at the tearing CMOD are also small and range from -2.35% to 3.36%. The per- 
cent difference between all of the solutions and the test load at tearing are also tabulated, and all but 
int. 5 are within 5% of the test result. 
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Figure 33. Influence of LPPL material choice on results for (a) P versus CMOD , (b) J((j)= 1 7°) 
versus CMOD , and (c) / versus (p. 
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Interpolation result int. 3 was chosen as representative of an ‘average’ choice for material 
properties and was then compared to the overall RR family of results for all laboratories. Figure 34 
shows the comparison of the P versus CMOD test data with the authors’ original FEA, the interpo- 
lated solution int. 3, and the analysis results of other RR participants, labs 2-15. The interpolated 
solution falls within the family of the RR results. It is important to recall that the interpolated solu- 
tion is generic. It provides a result out to the average CMOD n of its 16 nearest neighbors, so the hnal 
CMOD of the int. 3 solution is not part of the prediction; rather, it is the P versus CMOD trace 
up through a specified CMOD value that is of interest. Figure 35 compares the /-integral values at 
0=17° versus the CMOD results, the interpolated solution, and all of the labs’ solutions. Again, the 
interpolated solution is in excellent agreement with the family of solutions represented in the RR. 
Figure 36 shows the same data with the trajectory lines removed for all but the interpolated solu- 
tion. The RR participants were asked to provide /-integral results at three force levels (200, 252, and 
289 kN) as shown by the discrete ‘+’-shaped data points. To illustrate the reduced error when estab- 
lishing /-integral values from CMOD as opposed to force, each of the RR participants’ solutions 
was interpolated at CMOD values equal to the experimental CMOD at 256 kN and also at 0.15 mm. 
The reduced scatter in the / solutions is evident as revealed by the solid circles in figure 36. More- 
over, it is clear from the trace of the interpolated solution that it passes cleanly through the family of 
RR results at these CMOD values, providing an answer of equivalent quality to that expected from 
a custom finite element assessment of the test. 
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Figure 34. Interpolation, FEM, and RR participant results compared 
to experimental force versus CMOD response. 
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Figure 35. Comparison of interpolation, FEM, and RR participant 
results for J (0= 17°) versus CMOD. 
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36. Additional comparison of interpolation, FEM, and RR 
participant results for ./(</)= 17°) versus CMOD. 


4. CONCLUSIONS 


This set of elastic-plastic surface crack solutions and the interpolation methodologies devel- 
oped to sew them together into a complete, functional space represents a significant bridge for the 
practicing engineer toward commonplace assessment of surface cracks by the /-integral. This is par- 
ticularly true for the case of experimental evaluation of surface crack performance. The evaluation 
of fracture toughness in the elastic-plastic regime continues to lag significantly behind linear-elastic 
methods in implementation experimentally, even for common 2D geometries such as the C(T) speci- 
men, despite the advantages in flexibility and breadth of information elastic-plastic methods reveal. 
Surface crack toughness testing can provide the most direct measure of material performance in 
structurally representative configurations; yet, the detailed working knowledge of finite element mod- 
eling required to properly assess a surface crack test in the elastic-plastic regime (as most are) has kept 
surface crack toughness testing reserved mainly as a domain for the specialist. The costs in modeling 
time and software infrastructure are largely prohibitive for most experimental labs. The interpolation 
methodology described herein and solution space provided in the appendices represent a new evolu- 
tional step in tools for the analyst and experimentalist alike. 

Significant effort has been expended to ensure the reliability of this new tool. The authors have 
verified the solutions and the interpolation method through several techniques. The linear-elastic 
/-integral solutions were shown to be very similar to the Newman-Raju solutions. Domain conver- 
gence for the elastic-plastic /-integral values at the final load step was demonstrated for the complete 
set of solutions at all crack front locations. Twenty-five benchmark FEMs were created purposefully 
exploiting gaps in the geometry and material solution matrices to test the effectiveness of the interpo- 
lation method. The interpolated solutions were able to predict the benchmark /-integral and reaction 
force solutions for a given CMOD value to within a few percent. The interpolation method was also 
used to estimate a solution for an interlaboratory round robin problem, and the interpolated solution 
had excellent agreement with the multiple laboratory results. The results of these studies indicate that 
the authors’ surface crack solution set and interpolation method are sufficiently accurate to provide 
easily obtained and useful solutions for the /-integral and P-CMOD response. 

There remain a number of refinements and extensions to this work that would be equally 
valuable. The most obvious extension is to develop a similar solution space for surface cracks in 
flat plates subjected to a bending stress boundary condition. This is feasible with much of the same 
framework used to create the tension case. Within the existing tension space there are opportunities to 
fill gaps in the geometric or material variables to improve interpolation errors or to implement more 
sophisticated interpolation schemes. This may be a warranted investigation prior to developing the 
bending solutions. The solution set provides an excellent opportunity to investigate general surface 
crack behavior such as refining surface crack deformation limits in LEFM and EPFM, developing 
a robust predictor of surface crack tearing location in 0, or extrapolating solutions to higher deforma- 
tion levels. Despite the many possibilities, the most hopeful outcome for the authors is adoption in 
the experimental community. Perhaps these solutions will eventually provide an impact similar to the 
equations of Newman and Raju. 
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APPENDIX A— WIDTH EFFECTS SUBSTUDY 


A substudy was conducted to investigate the effects of increasing WI2c on the /-integral solu- 
tions and the force versus CMOD response. The same model used in the normalization case study was 
utilized here with W/2c = 5, L/W= 2, a/B = 0.6, ale- 0.6, n = 6, cr,, = 1, v=0.3, and E/o v =300. Two 
additional FEMs were created and analyzed with WI2c= 10 and 20 with all other geometric propor- 
tions and properties the same. In addition, interpolated solutions were estimated for all three plate 
geometries. Figure 37 illustrates the normalized far-held stress, a n , versus CMOD n for all the FEMs 
and interpolated solutions. All of the o n versus CMOD n interpolated solutions are identical and 
match the WI2c = 5 FEM result, and all of the responses are essentially the same in the linear-elastic 
region. The wider FEM plate results require a slightly higher far-held stress to reach a given CMOD n 
value in the nonlinear region. Figure 38 shows the hgure 37 results with the stress axis converted to 
applied force, P. As expected, the interpolated result and the W/2c = 5 FEM result match exactly. 
The interpolated solutions of the wider plate results are a reasonable match in the linear regime, but 
the interpolated results slightly underpredict the force for a given CMOD n in the nonlinear region. 
Matching the P versus CMOD n response is not required to interpolate to a /-integral solution, but 
is important for a comparison to P versus CMOD experimental results. Figure 39 illustrates that the 
/ w (0= 90°) versus CMOD n results are almost identical for the FEA and interpolated solutions show- 
ing that the /-integral as a function of CMOD is essentially invariant to increasing plate width. The 
maximum CMOD n for the WI2c = 20 FEM is shown as a vertical dashed line in hgure 39, and hgure 
40 shows J n versus 0 results at this CMOD n level for all the FEA and interpolated solutions. Again, 
the J n versus 0 results are demonstrated to be insensitive to plate width for a given CMOD n level, and 
the interpolated solutions are an excellent match to the FEA solutions. 
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Figure 37. Comparison of the o n versus CMOD n results for geometries 
with differing 2d W. 



Figure 38. Comparison of the P versus CMOD n results for geometries 
with differing 2d W. 
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Figure 39. Comparison of the J n ((j)= 90°) versus CMOD n results for 
geometries with differing 2 c/W. 



Figure 40. Comparison of the J n versus (j) results for geometries 
with differing 2d W at the maximum CMOD value of 
the 2d W=20 FEM. 


APPENDIX B— COLLECTION OF BENCHMARK MODELS COMPARISON PLOTS 
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0 0.01 0.02 0.03 0.04 0 15 30 45 60 75 

CMOD n (|) (deg) 

Figure Bl. Comparison of benchmark model No. 1 FEM results and interpolated results for (a) P versus CMOD ; 

(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD r 
at <p= 30° and 90°, (e) J n versus P at <p= 30° and 90°, and (f) J n versus (j ) for CMOD n values A and B. 
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0 0.005 0.01 0.015 0 15 30 45 60 75 

CMOD n c|) (deg) 

Figure B2. Comparison of benchmark model No. 2 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD r 
at (j)= 30° and 90°, (e) J n versus P at (j)= 30° and 90°, and (f) J n versus (j ) for CMOD n values A and B.. 
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Figure B3. Comparison of benchmark model No. 3 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD r 
at (j)= 30° and 90°, (e) J n versus P at <p= 30° and 90°, and (f) J n versus <p for CMOD n values A and B. 
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Figure B4. Comparison of benchmark model No. 4 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD r 
at (j)= 30° and 90°, (e) J n versus P at <p = 30° and 90°, and (f) J n versus 0for CMOD n values A and B. 
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Figure B5. Comparison of benchmark model No. 5 FEM results and interpolated results for (a) P versus CMOD ; 

(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD r 
at (j)= 30° and 90°, (e) J n versus P at (j) = 30° and 90°, and (f) J n versus 0for CMOD n values A and B. 
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Figure B6. Comparison of benchmark model No. 6 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD } 
at (j)= 30° and 90°, (e) J n versus P at <p = 30° and 90°, and (f) J n versus <p for CMOD n values A and B. 
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Figure B8. Comparison of benchmark model No. 8 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD r 
at (j)= 30° and 90°, (e) J n versus P at </)= 30° and 90°, and (f) J n versus b for CMOD n values A and B.. 
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Figure B9. Comparison of benchmark model No. 9 FEM results and interpolated results for (a) P versus CMOD i 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD } 
at (j) = 30° and 90°, (e) J n versus P at <p= 30° and 90°, and (f) J n versus 0for CMOD n values A and B. 



66 


CMOD n x10 - 3 <|> (deg) 

Figure BIO. Comparison of benchmark model No. 10 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD ]V (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at (j) = 30° and 90°, and (f) J n versus <p for CMOD n values A and B. 
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Figure Bll. Comparison of benchmark model No. 1 1 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at <p = 30° and 90°, and (f) J n versus 0for CMOD n values A and B. 
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Figure B12. Comparison of benchmark model No. 12 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at <p= 30° and 90°, (e) J n versus P at <p = 30° and 90°, and (f) J n versus b for CMOD n values A and B. 
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Figure B13. Comparison of benchmark model No. 13 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at <p = 30° and 90°, and (f) J n versus (j ) for CMOD n values A and B. 
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Figure B14. Comparison of benchmark model No. 14 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at <p = 30° and 90°, and (f) J n versus b for CMOD n values A and B. 
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Figure B15. Comparison of benchmark model No. 15 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at <p= 30° and 90°, and (1) J n versus b for CMOD n values A and B. 
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Figure B16. Comparison of benchmark model No. 16 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at (j)= 30° and 90°, and (f) J n versus d for CMOD n values A and B. 
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Figure B17. Comparison of benchmark model No. 17 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at <p= 30° and 90°, (e) J n versus P at (/) = 30° and 90°, and (f) J n versus (j ) for CMOD n values A and B. 
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Figure B18. Comparison of benchmark model No. 18 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at <p = 30° and 90°, and (f) J n versus (j ) for CMOD n values A and B. 
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Figure B19. Comparison of benchmark model No. 19 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at <p= 30° and 90°, (e) J n versus P at (j) = 30° and 90°, and (f) J n versus (j ) for CMOD n values A and B. 
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Figure B20. Comparison of benchmark model No. 20 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at (j) = 30° and 90°, and (f) J n versus (j ) for CMOD n values A and B. 
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Figure B21. Comparison of benchmark model No. 21 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at (j) = 30° and 90°, and (f) J n versus <p for CMOD n values A and B. 
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Figure B22. Comparison of benchmark model No. 22 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at 0 = 30° and 90°, and (f) J n versus 0for CMOD n values A and B. 
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Figure B23. Comparison of benchmark model No. 23 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at (j)= 30° and 90°, (e) J n versus P at <p= 30° and 90°, and (f) J n versus b for CMOD n values A and B. 
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Figure B24. Comparison of benchmark model No. 24 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at <p= 30° and 90°, (e) J n versus P at <p= 30° and 90°, and (f) J n versus b for CMOD n values A and B. 
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Figure B25. Comparison of benchmark model No. 25 FEM results and interpolated results for (a) P versus CMOD 
(b) M versus CMOD n , (c) benchmark FEM/interpolated results versus CMOD n , (d) J n versus CMOD n 
at 0= 30° and 90°, (e) J n versus P at <p= 30° and 90°, and (1) J n versus b for CMOD n values A and B. 


APPENDIX C— FINITE ELEMENT MODELS SOLUTION DATABASE 


The 600 finite element solutions developed for this work have been compiled into an easily 
indexed database for distribution. Each independent FEM solution database follows exactly the 
same format and has the same number of rows and columns. Table 4 shows an example of the data- 
base format for one model. Each model is given a unique ‘index’ name made up of the a/B, ale, n, 
and Elc 7 ys values joined by underscores, e.g. ‘0.20_0.20_3_100.’ The geometry and material property 
values are then listed under their appropriate labels. Each FEM result is divided into 20 even incre- 
ments (steps) of CMOD n based on the maximum CMOD n in the specific solution, and the /-integral, 
reaction force, and CMOD results are interpolated to these 20 increments. Step numbers are listed in 
the far left column followed by the corresponding normalized far-field stress, a n , and CMOD n val- 
ues. The remainder of the columns list the normalized /-integral values, J n , for 0 > 0 > 90 in 2° incre- 
ments of </> for each of the 20 load steps. The database for each individual FEM terminates with the 
keyword ‘end.’ All of the FEM solutions are collected into a single tab delimited text file that can be 
easily processed by computer programs. The database file is available upon request from the NASA 
Center for Aerospace Information (CASI) at <http://www.sti.nasa.gov>. 


Table 4. Example database format for FEM solutions 


Index 

a/B 

ale 

n 

E/Sys 

W 

B 

Sys 






0.20 0.20_3_1 00 

0.2 

0.2 

3 

100 

10 

1 

1 






** 

** 

phi 

0 

2 

4 

6 

8 

10 

12 


88 

90 

step 

stress 

CMOD 

JO 

J2 

J4 

J6 

J8 

J10 

J12 


J88 

J90 

1 

7.937E-02 

8.397E-04 

1.240E-05 

1.235E-05 

1.224E-05 

1.255E-05 

1.313E-05 

1.392E-05 

1.492E-05 


4.948E-05 

4.949E-05 

2 

1.586E-01 

1.679E-03 

4.581 E-05 

4.606E-05 

4.562E-05 

4.674E-05 

4.889E-05 

5.182E-05 

5.554E-05 


1.841 E-04 

1.84 IE-04 

3 

2.376E-01 

2.519E-03 

9.888E-05 

1.013E-04 

1.004E-04 

1.028E-04 

1.074E-04 

1.138E-04 

1.219E-04 


4.034E-04 

4.035E-04 

4 

3.161E-01 

3.359E-03 

1.695E-04 

1.780E-04 

1.772E-04 

1.81 IE-04 

1.890E-04 

2.000E-04 

2.142E-04 


7.070E-04 

7.072E-04 

5 

3.940E-01 

4.198E-03 

2.565E-04 

2.772E-04 

2.790E-04 

2.850E-04 

2.969E-04 

3.138E-04 

3.356E-04 


1.103E-03 

1.104E-03 

6 

4.711E-01 

5.038E-03 

3.566E-04 

3.955E-04 

4.045E-04 

4.144E-04 

4.31 IE-04 

4.547E-04 

4.854E-04 


1.586E-03 

1.587E-03 

7 

5.472E-01 

5.878E-03 

4.681 E-04 

5.302E-04 

5.516E-04 

5.685E-04 

5.917E-04 

6.231 E-04 

6.640E-04 


2.152E-03 

2.153E-03 

8 

6.222E-01 

6.717E-03 

5.907E-04 

6.798E-04 

7.185E-04 

7.467E-04 

7.787E-04 

8.199E-04 

8.724E-04 


2.801 E-03 

2.802E-03 

9 

6.956E-01 

7.557E-03 

7.244E-04 

8.449E-04 

9.046E-04 

9.479E-04 

9.923E-04 

1.046E-03 

1.112E-03 


3.525E-03 

3.526E-03 

10 

7.669E-01 

8.397E-03 

8.697E-04 

1.026E-03 

1.11 0E-03 

1.171 E-03 

1.231 E-03 

1.300E-03 

1.383E-03 


4.306E-03 

4.308E-03 

11 

8.320E-01 

9.236E-03 

1.021E-03 

1.215E-03 

1.325E-03 

1.406E-03 

1.482E-03 

1.569E-03 

1.670E-03 


5.148E-03 

5.150E-03 

12 

8.887E-01 

1.008E-02 

1.185E-03 

1.423E-03 

1.561E-03 

1.663E-03 

1.758E-03 

1.863E-03 

1.985E-03 


6.002E-03 

6.004E-03 

13 

9.355E-01 

1.092E-02 

1.354E-03 

1.641E-03 

1.809E-03 

1.935E-03 

2.049E-03 

2.175E-03 

2.319E-03 


6.872E-03 

6.874E-03 

14 

9.701 E-01 

1.176E-02 

1.513E-03 

1.848E-03 

2.046E-03 

2.196E-03 

2.330E-03 

2.477E-03 

2.643E-03 


7.754E-03 

7.755E-03 

15 

9.914E-01 

1.260E-02 

1.658E-03 

2.040E-03 

2.268E-03 

2.440E-03 

2.595E-03 

2.762E-03 

2.950E-03 


8.646E-03 

8.648E-03 

16 

1.012E+00 

1.343E-02 

1.797E-03 

2.225E-03 

2.481 E-03 

2.675E-03 

2.850E-03 

3.037E-03 

3.246E-03 


9.544E-03 

9.546E-03 

17 

1.030E+00 

1.427E-02 

1.917E-03 

2.389E-03 

2.670E-03 

2.885E-03 

3.077E-03 

3.282E-03 

3.512E-03 


1.046E-02 

1.046E-02 

18 

1.048E+00 

1.51 IE-02 

2.035E-03 

2.553E-03 

2.861 E-03 

3.096E-03 

3.307E-03 

3.532E-03 

3.783E-03 


1.138E-02 

1.139E-02 

19 

1.067E+00 

1.595E-02 

2.154E-03 

2.720E-03 

3.057E-03 

3.313E-03 

3.542E-03 

3.788E-03 

4.060E-03 


1.233E-02 

1.233E-02 

20 

1.086E+00 

1.679E-02 

2.274E-03 

2.889E-03 

3.255E-03 

3.532E-03 

3.781 E-03 

4.047E-03 

4.342E-03 


1.327E-02 

1.327E-02 

end 
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APPENDIX D— BENCHMARK FINITE ELEMENT MODELS SOLUTION DATABASE 


The 25 benchmark finite element solutions developed for this work have been complied into 
an easily indexed database for distribution. The solution database follows the same format as the 
600 model FEM database described in appendix C. The database hie is available upon request from 
the NASA Center for Aerospace Information (CASI) at <http://www.sti.nasa.gov>. 
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